<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Uniform sampling of E. coli core</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin-bottom: 20px; padding-bottom: 4px;  }
.S3 { margin: 0px; padding: 10px 0px 10px 5px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: start;  }
.S4 { margin: -1px 0px 0px; padding: 10px 0px 10px 7px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: start;  }
.S5 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S6 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S7 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S8 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }
.S9 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S16 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S17 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S18 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0' id = 'T_D3C2B58F' ><span>Uniform sampling of E. coli core</span></h1><h2  class = 'S1' id = 'H_1BBB3AA5' ><span>Author(s): German A. Preciat Gonzalez, Ronan M.T. Fleming, Leiden University.</span><span style=' font-weight: bold;'>  </span></h2><h2  class = 'S1' id = 'H_D9B46067' ><span>Reviewer(s): </span></h2><div  class = 'S2'><div  class = 'S3'><span style=' font-weight: bold;'>Table of Contents</span></div><div  class = 'S4'><a href = "#H_1BBB3AA5"><span>Author(s): German A. Preciat Gonzalez, Ronan M.T. Fleming, Leiden University.  
</span></a><a href = "#H_D9B46067"><span>Reviewer(s): 
</span></a><a href = "#H_8DFA11B4"><span>INTRODUCTION
</span></a><a href = "#H_7E2A567B"><span>MATERIALS - EQUIPMENT SETUP
</span></a><a href = "#H_B642E8E4"><span>PROCEDURE
</span></a><a href = "#H_ED106D18"><span>Load E. coli core model
</span></a><a href = "#H_2E9842A7"><span>E. coli core in aerobic an anaerobic conditions
</span></a><a href = "#H_F31F83BB"><span>Flux variability analysis
</span></a><a href = "#H_FCF7B08E"><span>Undersampling
</span></a><a href = "#H_4781AD3B"><span>Sampling
</span></a><a href = "#H_D53645CF"><span>Acknowledgements
</span></a><a href = "#H_67123C98"><span>References</span></a></div></div><h2  class = 'S1' id = 'H_8DFA11B4' ><span>INTRODUCTION</span></h2><div  class = 'S5'><span>The flux space </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: normal; font-weight: normal; color: rgb(0, 0, 0);">Ω</span><span> for a given set of biochemical and physiologic constraints is represented by: </span></div><div  class = 'S6'><span texencoding="\Omega = \{v \mid Sv=b; l \leq v\leq u\}" style="vertical-align:-5px"><img src="" width="168.5" height="19" /></span></div><div  class = 'S5'><span>where </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">v</span><span> represents feasible flux vectors, </span><span> </span><span texencoding="S\in\mathcal{Z}^{m\times n}" style="vertical-align:-5px"><img src="" width="63.5" height="19" /></span><span> the stoichiometric matrix, while </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">l</span><span> and </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">u</span><span> are lower and upper bounds on fluxes. These criteria still allow a wide range of admissible flux distributions which, in FBA are commonly further restricted by introducing an objective to optimise, transforming the question of admissible fluxes into an FBA problem</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span> of the form</span></div><div  class = 'S6'><span texencoding="\begin{array}{ll}
\min\limits _{v} &amp; c^{T}v\\
\text{s.t.} &amp; Sv=b,\\
 &amp; l\leq v\leq u,
\end{array}" style="vertical-align:-30px"><img src="" width="97.5" height="72" /></span></div><div  class = 'S5'><span>where </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">c</span><span> is a linear biological objective function (biomass, ATP consumption, HEME production, etc.). Even under these conditions there is commonly a range of optimal flux distributions, which can be investigated using flux variability analysis. If the general capabilities of the model are of interest, however, uniform sampling of the entire flux space</span><span> </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: normal; font-weight: normal; color: rgb(0, 0, 0);">Ω</span><span> is able to provide an unbiased characterization, and therefore, can be used to investigate the biochemical networks. It requires collecting a statistically meaningful number of flux distributions uniformly spread throughout the whole flux space and then analysing their properties. There are three basic steps to perform a uniform sampling for a set of feasible fluxes:</span></div><ul  class = 'S7'><li  class = 'S8'><span>Define the flux space to be sampled from physical and biochemical constraints</span></li><li  class = 'S8'><span>Randomly sample the defined flux space based on uniform statistical criteria</span></li><li  class = 'S8'><span>If is necessary, section the flux space according to post-sampling.</span></li></ul><div  class = 'S5'><span>In COBRA v3 the default sampling algorithm is coordinate hit-and-run with rounding (CHRR). This algorithm first rounds the anisotropic flux space </span><span> </span><span>Ω</span><span> using a maximum volume ellipsoid algorithm</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;3&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span> and then performs a uniform sampling based on the provably efficient hit-and-run random walk</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;4&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>. Below is a high-level illustration of the process to uniformly sample a random metabolic flux vector </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">v</span><span> from the set </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: normal; font-weight: normal; color: rgb(0, 0, 0);">Ω</span><span> of all feasible metabolic fluxes (grey). </span><span style=' font-weight: bold;'>1)</span><span> Apply a rounding transformation </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">T</span><span> to </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: normal; font-weight: normal; color: rgb(0, 0, 0);">Ω</span><span>. The transformed set </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi&gt;&amp;ohm;&lt;/mi&gt;&lt;mo&gt;&amp;prime;&lt;/mo&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;T&lt;/mi&gt;&lt;mi&gt;&amp;ohm;&lt;/mi&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="61" height="18" /></span><span> is such that its maximal inscribed ellipsoid (blue) approximates a unit ball. </span><span style=' font-weight: bold;'>2)</span><span> Take </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">q</span><span> steps of coordinate hit-and-run. At each step, i) pick a random coordinate direction </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;e&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;i&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="12.5" height="20" /></span><span>, and ii) move from current point </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;v&lt;/mi&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mo&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mo&gt;&amp;isinv;&lt;/mo&gt;&lt;mi&gt;&amp;ohm;&lt;/mi&gt;&lt;mo&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="53" height="20" /></span><span> to a random point </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;v&lt;/mi&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mo stretchy=&quot;false&quot;&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;mo&gt;+&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mo stretchy=&quot;false&quot;&gt;&amp;isinv;&lt;/mo&gt;&lt;mi&gt;&amp;ohm;&lt;/mi&gt;&lt;mo stretchy=&quot;false&quot;&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="65.5" height="20" /></span><span> along </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;v&lt;/mi&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mo&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mo&gt;+&lt;/mo&gt;&lt;mi&gt;&amp;alpha;&lt;/mi&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;e&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;i&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mo&gt;&amp;cap;&lt;/mo&gt;&lt;mi&gt;&amp;ohm;&lt;/mi&gt;&lt;mo&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="85" height="20" /></span><span>. </span><span style=' font-weight: bold;'>3)</span><span> Map samples back to the original space by applying the inverse</span><span> transformation, i.e., </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;v&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;T&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mo&gt;&amp;minus;&lt;/mo&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;v&lt;/mi&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mo&gt;&amp;prime;&lt;/mo&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;k&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="70.5" height="21" /></span><span>.</span></div><div  class = 'S6'><img class = "imageNode" src = "" width = "647" height = "280" alt = "" style = "vertical-align: baseline"></img></div><h2  class = 'S9' id = 'H_7E2A567B' ><span>MATERIALS - EQUIPMENT SETUP</span></h2><div  class = 'S5'><span>Please ensure that all the required dependencies (e.g. , </span><span style=' font-family: monospace;'>git</span><span> and </span><span style=' font-family: monospace;'>curl</span><span>) of The COBRA Toolbox have been properly installed by following the installation guide </span><a href = "https://opencobra.github.io/cobratoolbox/stable/installation.html"><span>here</span></a><span>. Please ensure that the COBRA Toolbox has been initialised (tutorial_initialize.mlx) and verify that the pre-packaged LP and QP solvers are functional (tutorial_verify.mlx).</span></div><div  class = 'S5'><span>Please note that some of the plotting options in the tutorial require Matlab 2016a or higher. Moreover, the tutorial requires a working installation of the Parallel Computing Toolbox.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% uncomment this line below to see what toolboxes are installed</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% ver</span></span></div></div></div><div  class = 'S12'><span>In this tutorial, we will perform FVA using the function </span><span style=' font-family: monospace;'>fluxVariability</span><span>. Change the variable </span><span style=' font-family: monospace;'>options.useFastFVA = 1</span><span> to use </span><span style=' font-family: monospace;'>fastFVA </span><span>instead</span><span style=' font-family: monospace;'>.</span><span> Note that the solver ibm_cplex is required for the function fastFVA.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S13' id = 'H_636DC3ED' ><span style="white-space: pre"><span >options.useFastFVA = 0;</span></span></div></div></div><h2  class = 'S1' id = 'H_B642E8E4' ><span>PROCEDURE</span></h2><h2  class = 'S1' id = 'H_ED106D18' ><span>Load E. coli core model</span></h2><div  class = 'S5'><span>The most appropriate way to load a model into The COBRA Toolbox is to use the </span><span style=' font-family: monospace;'>readCbModel</span><span> function. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >fileName = </span><span style="color: rgb(170, 4, 249);">'ecoli_core_model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >~exist(</span><span style="color: rgb(170, 4, 249);">'modelOri'</span><span >,</span><span style="color: rgb(170, 4, 249);">'var'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    modelOri = readCbModel(fileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%backward compatibility with primer requires relaxation of upper bound on</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%ATPM</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >modelOri = changeRxnBounds(modelOri,</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >,1000,</span><span style="color: rgb(170, 4, 249);">'u'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >model = modelOri;</span></span></div></div></div><h2  class = 'S1' id = 'H_2E9842A7' ><span>E. coli core in aerobic an anaerobic conditions</span></h2><div  class = 'S5'><span>Remove the objective from the model and set a small lower bound on the rate of the biomass reaction</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >biomassRxnAbbr = </span><span style="color: rgb(170, 4, 249);">'Biomass_Ecoli_core_N(w/GAM)-Nmet2'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >ibm = find(ismember(model.rxns, biomassRxnAbbr));  </span><span style="color: rgb(2, 128, 9);">% column index of the biomass reaction</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >model.lb(ibm)=0.05;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >model.c(:)=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'></div></div></div><div  class = 'S12'><span>We will investigate ATP energy production with limited and unlimited oxygen uptake. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >aerobicModel = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >,-17,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >anAerobicModel = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >,-1,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><h2  class = 'S9' id = 'H_F31F83BB' ><span>Flux variability analysis</span></h2><div  class = 'S5'><span>Flux variability analysis (FVA) returns the minimum and maximum possible flux through every reaction in a model.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >options.useFastFVA</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    [minUn, maxUn] = fastFVA(aerobicModel, 100);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    [minLim, maxLim] = fastFVA(anAerobicModel, 100);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    [minUn, maxUn] = fluxVariability(aerobicModel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    [minLim, maxLim] = fluxVariability(anAerobicModel);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S12'><span></span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'Max. biomass production with oxygen uptake: %.4f/h.\n'</span><span >, maxUn(ibm));</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="ADBD9BF4" data-testid="output_0" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Max. biomass production with oxygen uptake: 0.7637/h.</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S17'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'Max. biomass production without oxygen uptake: %.4f/h.\n\n'</span><span >, maxLim(ibm));</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="089CDE9E" data-testid="output_1" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Max. biomass production without oxygen uptake: 0.2477/h.</div></div></div></div></div><div  class = 'S12'><span>An overall comparison of the FVA results can be obtained by computing the </span><a href = "https://en.wikipedia.org/wiki/Jaccard_index"><span>Jaccard index</span></a><span> for each reaction. The Jaccard index is here defined as the ratio between the intersection and union of the flux ranges in the aerobic and anaerobic models. A Jaccard index of 0 indicates completely disjoint flux ranges and a Jaccard index of 1 indicates completely overlapping flux ranges. The mean Jaccard index gives an indication of the overall similarity between the models.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >J = fvaJaccardIndex([minUn, minLim], [maxUn, maxLim]);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S18'><span style="white-space: pre"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'Mean Jaccard index = %.4f.\n'</span><span >, mean(J));</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="38B3A00A" data-testid="output_2" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Mean Jaccard index = 0.5524.</div></div></div></div></div><div  class = 'S12'><span>To visualise the FVA results, we plot the flux ranges as errorbars, with reactions sorted by the Jaccard index.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >E = [(maxUn - minUn)/2 (maxLim - minLim)/2];</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >Y = [minUn minLim] + E;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >X = [(1:length(Y)) - 0.1; (1:length(Y)) + 0.1]';</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >[~, xj] = sort(J);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >f1 = figure;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >strcmp(version(</span><span style="color: rgb(170, 4, 249);">'-release'</span><span >), </span><span style="color: rgb(170, 4, 249);">'2016b'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    errorbar(X, Y(xj, :), E(xj, :), </span><span style="color: rgb(170, 4, 249);">'linestyle'</span><span >, </span><span style="color: rgb(170, 4, 249);">'none'</span><span >, </span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2, </span><span style="color: rgb(170, 4, 249);">'capsize'</span><span >, 0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    </span><span style="color: rgb(2, 128, 9);">%errorbar(X, Y(xj, :), E(xj, :), 'linestyle', 'none', 'linewidth', 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    errorbar(X(:,1), Y(xj, 1), E(xj, 1), </span><span style="color: rgb(170, 4, 249);">'linestyle'</span><span >, </span><span style="color: rgb(170, 4, 249);">'none'</span><span >, </span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2,</span><span style="color: rgb(170, 4, 249);">'Color'</span><span >,</span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >    errorbar(X(:,2), Y(xj, 2), E(xj, 2), </span><span style="color: rgb(170, 4, 249);">'linestyle'</span><span >, </span><span style="color: rgb(170, 4, 249);">'none'</span><span >, </span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2,</span><span style="color: rgb(170, 4, 249);">'Color'</span><span >,</span><span style="color: rgb(170, 4, 249);">'r'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >set(gca, </span><span style="color: rgb(170, 4, 249);">'xlim'</span><span >, [0, length(Y) + 1])</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Reaction'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Flux range (mmol/gDW/h)'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >ylim([-50,50])</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >yyaxis </span><span style="color: rgb(170, 4, 249);">right</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >plot(J(xj),</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >legend(</span><span style="color: rgb(170, 4, 249);">'Aerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Anaerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Jaccard'</span><span >,</span><span style="color: rgb(170, 4, 249);">'location'</span><span >, </span><span style="color: rgb(170, 4, 249);">'northoutside'</span><span >, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >       </span><span style="color: rgb(170, 4, 249);">'orientation'</span><span >, </span><span style="color: rgb(170, 4, 249);">'horizontal'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S18'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Jaccard index'</span><span >)</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="78F8C1F4" data-testid="output_3" style="width: 458px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><h2  class = 'S9' id = 'H_FCF7B08E' ><span>Undersampling</span></h2><div  class = 'S5'><span>CHRR can be called via the function </span><span style=' font-family: monospace;'>sampleCbModel</span><span>. The main inputs to </span><span style=' font-family: monospace;'>sampleCbModel</span><span> are a COBRA model structure, the name of the selected sampler and a parameter struct that controls properties of the sampler used. In the instance of CHRR, two parameters are important: the sampling density (</span><span style=' font-family: monospace;'>nStepsPerPoint)</span><span> and the number of samples (</span><span style=' font-family: monospace;'>nPointsReturned). </span><span>The total length of the random walk is </span><span style=' font-family: monospace;'>nStepsPerPoint*nPointsReturned</span><span>. The time it takes to run the sampler depends on the total length of the random walk and the size of the model</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>. However, using sampling parameters that are too small will lead to invalid sampling distributions, e.g.,</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.nStepsPerPoint = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >options.nPointsReturned = 500;</span></span></div></div></div><div  class = 'S12'><span>An additional on/off parameter (</span><span style=' font-family: monospace;'>toRound</span><span>) controls whether or not the polytope is rounded. Rounding large models can be slow but is strongly recommended for the first round of sampling. Below we show how to get around this step in subsequent rounds. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: pre"><span >options.toRound = 1;</span></span></div></div></div><div  class = 'S12'><span>The method outputs two results. First, the model used for sampling (in case of </span><span style=' font-family: monospace;'>toRound = 1</span><span> this would be the rounded model), and second, the samples generated. </span></div><div  class = 'S5'><span>To sample the aerobic and anaerobic E. coli core models use:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: pre"><span >[P_un, X1_un] =  sampleCbModel(aerobicModel, [], [], options);</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="F6164B2C" data-testid="output_4" data-width="428" data-height="59" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Checking for width 0 facets...
Currently (P.A, P.b) are in 95 dimensions
Checking for width 0 facets...
Found 8 degenerate reactions, adding them to the equality subspace.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="5A9A3101" data-testid="output_5" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Rank deficient, rank = 71, tol =  1.052324e-10.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="3379EEAF" data-testid="output_6" data-width="428" data-height="101" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Now in 24 dimensions after restricting
Removed 174 zero rows
Preconditioning A with gmscale
Rounding...
  Converged!
Maximum volume ellipsoid found, and the origin is inside the transformed polytope.
Generating samples...</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S17'><span style="white-space: pre"><span >[P_lim, X1_lim] = sampleCbModel(anAerobicModel, [], [], options);</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="05AA6E0C" data-testid="output_7" data-width="428" data-height="59" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Checking for width 0 facets...
Currently (P.A, P.b) are in 95 dimensions
Checking for width 0 facets...
Found 8 degenerate reactions, adding them to the equality subspace.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsWarningElement" uid="C8116ACB" data-testid="output_8" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType" style="white-space: normal; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"><div class="diagnosticMessage-messagePart" style="white-space: pre-wrap; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;">Warning: Rank deficient, rank = 71, tol =  1.052324e-10.</div><div class="diagnosticMessage-stackPart" style="white-space: pre; font-style: normal; color: rgb(255, 100, 0); font-size: 12px;"></div></div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="BA57D644" data-testid="output_9" data-width="428" data-height="101" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Now in 24 dimensions after restricting
Removed 174 zero rows
Preconditioning A with gmscale
Rounding...
  Converged!
Maximum volume ellipsoid found, and the origin is inside the transformed polytope.
Generating samples...</div></div></div></div></div><div  class = 'S5'><span>The sampler outputs the sampled flux distributions (X_un and X_lim) and the rounded polytope (P_un and P_lim). Histograms of sampled ATP synthase show that the models are severely undersampled, as evidenced by the presence of multiple sharp peaks.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >nbins = 20;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >[yUn, xUn] = hist(X1_un(ibm, :), nbins,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >[yLims, xLims] = hist(X1_lim(ibm, :), nbins,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >f2 = figure;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >plot(xUn, yUn, xLims, yLims,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >legend(</span><span style="color: rgb(170, 4, 249);">'Aerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Anaerobic'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Biomass Flux (mmol/gDW/h)'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S18'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'# samples'</span><span >)</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="6B5727B5" data-testid="output_10" style="width: 458px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S5'><span style=' font-weight: bold;'>Figure 1. </span><span>Undersampling results from selecting too small sampling parameters. </span></div><h2  class = 'S1' id = 'H_4781AD3B' ><span>Sampling</span></h2><div  class = 'S5'><span>The appropriate parameter values depend on the dimension of the polytope Ω  defined by the model constraints (see intro). One rule of thumb says to set  </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;nSkip&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mn&gt;8&lt;/mn&gt;&lt;mo&gt;*&lt;/mo&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;dim&lt;/mi&gt;&lt;mrow&gt;&lt;mo&gt;(&lt;/mo&gt;&lt;mrow&gt;&lt;mi&gt;&amp;ohm;&lt;/mi&gt;&lt;/mrow&gt;&lt;mo&gt;)&lt;/mo&gt;&lt;/mrow&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="130" height="20" /></span><span> to ensure the statistical independence of samples. The random walk should be long enough to ensure convergence to a stationary sampling distribution</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.nStepsPerPoint = 8 * size(P_lim.A, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span >options.nPointsReturned = 1000;</span></span></div></div></div><div  class = 'S12'><span>This time, we can avoid the rounding step by inputting the rounded polytope from the previous round of sampling.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >options.toRound = 0;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S18'><span style="white-space: pre"><span >[~, X2_un] = sampleCbModel(aerobicModel, [], [], options, P_un);</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="02BC0643" data-testid="output_11" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Generating samples...</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S17'><span style="white-space: pre"><span >[~, X2_lim] = sampleCbModel(anAerobicModel, [], [], options, P_lim);</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="2016836E" data-testid="output_12" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Generating samples...</div></div></div></div></div><div  class = 'S5'><span>The converged sampling distributions for the biomass reaction are much smoother, with a single peak near the minimum biomass flux.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >nbins = 20;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >[yUn, xUn] = hist(X2_un(ibm, :), nbins,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >[yLims, xLims] = hist(X2_lim(ibm, :), nbins,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >f3 = figure;</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >p1 = plot(xUn, yUn, xLims, yLims,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >legend(</span><span style="color: rgb(170, 4, 249);">'Aerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Anaerobic'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Biomass Flux (mmol/gDW/h)'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S18'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'# samples'</span><span >)</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="4A4AE90B" data-testid="output_13" style="width: 458px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S12'><span>Adding </span><span>the FVA res</span><span>u</span><span>lts to the plot shows that the sampling distributions give more detailed information about the differences between the two models. In particular, we see that the biomass flux minima and maxima are not equally probable. The number of samples from both the aerobic and anaerobic models peaks at the minimum flux of zero, and decreases approximately monotonically towards the maximum growth rate. It decreases more slowly in the unanaerobic model, indicating that higher ATP production is more probable under unlimited oxygen uptake conditions. It is interesting to see that maximum ATP production is highly improbable in both models.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >ylims = get(gca, </span><span style="color: rgb(170, 4, 249);">'ylim'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >cUn = get(p1(1), </span><span style="color: rgb(170, 4, 249);">'color'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >cLim = get(p1(2), </span><span style="color: rgb(170, 4, 249);">'color'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">on</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >p2 = plot([minUn(ibm), minUn(ibm)], ylims, </span><span style="color: rgb(170, 4, 249);">'--'</span><span >, [maxUn(ibm), maxUn(ibm)], ylims, </span><span style="color: rgb(170, 4, 249);">'--'</span><span >,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >set(p2,</span><span style="color: rgb(170, 4, 249);">'color'</span><span >, cUn)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >p3 = plot([minLim(ibm), minLim(ibm)], ylims, </span><span style="color: rgb(170, 4, 249);">'--'</span><span >, [maxLim(ibm), maxLim(ibm)], ylims, </span><span style="color: rgb(170, 4, 249);">'--'</span><span >,</span><span style="color: rgb(170, 4, 249);">'linewidth'</span><span >, 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >set(p3, </span><span style="color: rgb(170, 4, 249);">'color'</span><span >, cLim)</span></span></div></div><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: pre"><span >legend(</span><span style="color: rgb(170, 4, 249);">'Aerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'Anaerobic'</span><span >, </span><span style="color: rgb(170, 4, 249);">'AerobicFVA'</span><span >,</span><span style="color: rgb(170, 4, 249);">'AerobicFVA'</span><span >, </span><span style="color: rgb(170, 4, 249);">'AnaerobicFVA'</span><span >,</span><span style="color: rgb(170, 4, 249);">'AnaerobicFVA'</span><span >)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S18'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">off</span></span></div><div  class = 'S16'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="0D14D853" data-testid="output_14" style="width: 458px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S12'><span></span></div><h2  class = 'S9' id = 'H_D53645CF' ><span>Acknowledgements</span></h2><div  class = 'S5'><span>Based on a genome-scale sampling tutorial by Hulda S. Haraldsdóttir and German A. Preciat Gonzalez.</span></div><h2  class = 'S1' id = 'H_67123C98' ><span>References</span></h2><div  class = 'S5'><span>1. Orth, J. D., Thiele I., and Palsson, B. Ø. What is flux balance analysis? </span><span style=' font-style: italic;'>Nat. Biotechnol.</span><span> 28(3), 245-248 (2010).</span></div><div  class = 'S5'><span>2. Haraldsdóttir, H. S., Cousins, B., Thiele, I., Fleming, R.M.T., and Vempala, S. CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based metabolic models. </span><span style=' font-style: italic;'>Bioinformatics</span><span>. 33(11), 1741-1743 (2016).</span></div><div  class = 'S5'><span>3. Zhang, Y. and Gao, L. On Numerical Solution of the Maximum Volume Ellipsoid Problem. </span><span style=' font-style: italic;'>SIAM J. Optimiz</span><span>. 14(1), 53-76 (2001).</span></div><div  class = 'S5'><span>4. Berbee, H. C. P., Boender, C. G. E., Rinnooy Ran, A. H. G., Scheffer, C. L., Smith, R. L., Telgen, J. Hit-and-run algorithms for the identification of nonredundant linear inequalities. </span><span style=' font-style: italic;'>Math. Programming</span><span>, 37(2), 184-207 (1987).</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Uniform sampling of E. coli core
%% Author(s): German A. Preciat Gonzalez, Ronan M.T. Fleming, Leiden University.  
%% Reviewer(s): 
%% INTRODUCTION
% The flux space $\Omega$ for a given set of biochemical and physiologic constraints 
% is represented by: 
% 
% $$\Omega = \{v \mid Sv=b; l \leq v\leq u\}$$
% 
% where $v$ represents feasible flux vectors,  $S\in\mathcal{Z}^{m\times n}$ 
% the stoichiometric matrix, while $l$ and $u$ are lower and upper bounds on fluxes. 
% These criteria still allow a wide range of admissible flux distributions which, 
% in FBA are commonly further restricted by introducing an objective to optimise, 
% transforming the question of admissible fluxes into an FBA problem${\;}^1$ of 
% the form
% 
% $$\begin{array}{ll}\min\limits _{v} & c^{T}v\\\text{s.t.} & Sv=b,\\ & l\leq 
% v\leq u,\end{array}$$
% 
% where $c$ is a linear biological objective function (biomass, ATP consumption, 
% HEME production, etc.). Even under these conditions there is commonly a range 
% of optimal flux distributions, which can be investigated using flux variability 
% analysis. If the general capabilities of the model are of interest, however, 
% uniform sampling of the entire flux space $\Omega$ is able to provide an unbiased 
% characterization, and therefore, can be used to investigate the biochemical 
% networks. It requires collecting a statistically meaningful number of flux distributions 
% uniformly spread throughout the whole flux space and then analysing their properties. 
% There are three basic steps to perform a uniform sampling for a set of feasible 
% fluxes:
%% 
% * Define the flux space to be sampled from physical and biochemical constraints
% * Randomly sample the defined flux space based on uniform statistical criteria
% * If is necessary, section the flux space according to post-sampling.
%% 
% In COBRA v3 the default sampling algorithm is coordinate hit-and-run with 
% rounding (CHRR). This algorithm first rounds the anisotropic flux space  Ω using 
% a maximum volume ellipsoid algorithm${\;}^3$ and then performs a uniform sampling 
% based on the provably efficient hit-and-run random walk${\;}^4$. Below is a 
% high-level illustration of the process to uniformly sample a random metabolic 
% flux vector $v$ from the set $\Omega$ of all feasible metabolic fluxes (grey). 
% *1)* Apply a rounding transformation $T$ to $\Omega$. The transformed set $\Omega 
% \prime =T\Omega \;$ is such that its maximal inscribed ellipsoid (blue) approximates 
% a unit ball. *2)* Take $q$ steps of coordinate hit-and-run. At each step, i) 
% pick a random coordinate direction $e_i$, and ii) move from current point $v\prime_k 
% \in \Omega \prime$ to a random point $v\prime_{k+1} \in \Omega \prime$ along 
% $v\prime_k +\alpha e_i \cap \Omega \prime$. *3)* Map samples back to the original 
% space by applying the inverse transformation, i.e., $v_k =T^{-1} v\prime_k$.
% 
% 
%% MATERIALS - EQUIPMENT SETUP
% Please ensure that all the required dependencies (e.g. , |git| and |curl|) 
% of The COBRA Toolbox have been properly installed by following the installation 
% guide <https://opencobra.github.io/cobratoolbox/stable/installation.html here>. 
% Please ensure that the COBRA Toolbox has been initialised (tutorial_initialize.mlx) 
% and verify that the pre-packaged LP and QP solvers are functional (tutorial_verify.mlx).
% 
% Please note that some of the plotting options in the tutorial require Matlab 
% 2016a or higher. Moreover, the tutorial requires a working installation of the 
% Parallel Computing Toolbox.

% uncomment this line below to see what toolboxes are installed
% ver
%% 
% In this tutorial, we will perform FVA using the function |fluxVariability|. 
% Change the variable |options.useFastFVA = 1| to use |fastFVA| instead|.| Note 
% that the solver ibm_cplex is required for the function fastFVA.

options.useFastFVA = 0;
%% PROCEDURE
%% Load E. coli core model
% The most appropriate way to load a model into The COBRA Toolbox is to use 
% the |readCbModel| function. 

fileName = 'ecoli_core_model.mat';
if ~exist('modelOri','var')
    modelOri = readCbModel(fileName);
end
%backward compatibility with primer requires relaxation of upper bound on
%ATPM
modelOri = changeRxnBounds(modelOri,'ATPM',1000,'u');
model = modelOri;
%% E. coli core in aerobic an anaerobic conditions
% Remove the objective from the model and set a small lower bound on the rate 
% of the biomass reaction

biomassRxnAbbr = 'Biomass_Ecoli_core_N(w/GAM)-Nmet2';
ibm = find(ismember(model.rxns, biomassRxnAbbr));  % column index of the biomass reaction
model.lb(ibm)=0.05;
model.c(:)=0;

%% 
% We will investigate ATP energy production with limited and unlimited oxygen 
% uptake. 

aerobicModel = changeRxnBounds(model,'EX_o2(e)',-17,'l');
anAerobicModel = changeRxnBounds(model,'EX_o2(e)',-1,'l');
%% Flux variability analysis
% Flux variability analysis (FVA) returns the minimum and maximum possible flux 
% through every reaction in a model.

if options.useFastFVA
    [minUn, maxUn] = fastFVA(aerobicModel, 100);
    [minLim, maxLim] = fastFVA(anAerobicModel, 100);
else
    [minUn, maxUn] = fluxVariability(aerobicModel);
    [minLim, maxLim] = fluxVariability(anAerobicModel);
end
%% 
% 

fprintf('Max. biomass production with oxygen uptake: %.4f/h.\n', maxUn(ibm));
fprintf('Max. biomass production without oxygen uptake: %.4f/h.\n\n', maxLim(ibm));
%% 
% An overall comparison of the FVA results can be obtained by computing the 
% <https://en.wikipedia.org/wiki/Jaccard_index Jaccard index> for each reaction. 
% The Jaccard index is here defined as the ratio between the intersection and 
% union of the flux ranges in the aerobic and anaerobic models. A Jaccard index 
% of 0 indicates completely disjoint flux ranges and a Jaccard index of 1 indicates 
% completely overlapping flux ranges. The mean Jaccard index gives an indication 
% of the overall similarity between the models.

J = fvaJaccardIndex([minUn, minLim], [maxUn, maxLim]);
fprintf('Mean Jaccard index = %.4f.\n', mean(J));
%% 
% To visualise the FVA results, we plot the flux ranges as errorbars, with reactions 
% sorted by the Jaccard index.

E = [(maxUn - minUn)/2 (maxLim - minLim)/2];
Y = [minUn minLim] + E;
X = [(1:length(Y)) - 0.1; (1:length(Y)) + 0.1]';

[~, xj] = sort(J);

f1 = figure;
if strcmp(version('-release'), '2016b')
    errorbar(X, Y(xj, :), E(xj, :), 'linestyle', 'none', 'linewidth', 2, 'capsize', 0);
else
    %errorbar(X, Y(xj, :), E(xj, :), 'linestyle', 'none', 'linewidth', 2);
    hold on
    errorbar(X(:,1), Y(xj, 1), E(xj, 1), 'linestyle', 'none', 'linewidth', 2,'Color','b');
    errorbar(X(:,2), Y(xj, 2), E(xj, 2), 'linestyle', 'none', 'linewidth', 2,'Color','r');
end
set(gca, 'xlim', [0, length(Y) + 1])

xlabel('Reaction')
ylabel('Flux range (mmol/gDW/h)')
ylim([-50,50])
yyaxis right
plot(J(xj),'linewidth', 2)
legend('Aerobic', 'Anaerobic', 'Jaccard','location', 'northoutside', ...
       'orientation', 'horizontal')
ylabel('Jaccard index')
%% Undersampling
% CHRR can be called via the function |sampleCbModel|. The main inputs to |sampleCbModel| 
% are a COBRA model structure, the name of the selected sampler and a parameter 
% struct that controls properties of the sampler used. In the instance of CHRR, 
% two parameters are important: the sampling density (|nStepsPerPoint)| and the 
% number of samples (|nPointsReturned).| The total length of the random walk is 
% |nStepsPerPoint*nPointsReturned|. The time it takes to run the sampler depends 
% on the total length of the random walk and the size of the model${\;}^2$. However, 
% using sampling parameters that are too small will lead to invalid sampling distributions, 
% e.g.,

options.nStepsPerPoint = 1;
options.nPointsReturned = 500;
%% 
% An additional on/off parameter (|toRound|) controls whether or not the polytope 
% is rounded. Rounding large models can be slow but is strongly recommended for 
% the first round of sampling. Below we show how to get around this step in subsequent 
% rounds. 

options.toRound = 1;
%% 
% The method outputs two results. First, the model used for sampling (in case 
% of |toRound = 1| this would be the rounded model), and second, the samples generated. 
% 
% To sample the aerobic and anaerobic E. coli core models use:

[P_un, X1_un] =  sampleCbModel(aerobicModel, [], [], options);
[P_lim, X1_lim] = sampleCbModel(anAerobicModel, [], [], options);
%% 
% The sampler outputs the sampled flux distributions (X_un and X_lim) and the 
% rounded polytope (P_un and P_lim). Histograms of sampled ATP synthase show that 
% the models are severely undersampled, as evidenced by the presence of multiple 
% sharp peaks.

nbins = 20;
[yUn, xUn] = hist(X1_un(ibm, :), nbins,'linewidth', 2);
[yLims, xLims] = hist(X1_lim(ibm, :), nbins,'linewidth', 2);

f2 = figure;
plot(xUn, yUn, xLims, yLims,'linewidth', 2);
legend('Aerobic', 'Anaerobic')
xlabel('Biomass Flux (mmol/gDW/h)')
ylabel('# samples')
%% 
% *Figure 1.* Undersampling results from selecting too small sampling parameters. 
%% Sampling
% The appropriate parameter values depend on the dimension of the polytope Ω  
% defined by the model constraints (see intro). One rule of thumb says to set  
% $\textrm{nSkip}=8*{\dim \left(\Omega \right)}^2$ to ensure the statistical independence 
% of samples. The random walk should be long enough to ensure convergence to a 
% stationary sampling distribution${\;}^2$.

options.nStepsPerPoint = 8 * size(P_lim.A, 2);
options.nPointsReturned = 1000;
%% 
% This time, we can avoid the rounding step by inputting the rounded polytope 
% from the previous round of sampling.

options.toRound = 0;
[~, X2_un] = sampleCbModel(aerobicModel, [], [], options, P_un);
[~, X2_lim] = sampleCbModel(anAerobicModel, [], [], options, P_lim);
%% 
% The converged sampling distributions for the biomass reaction are much smoother, 
% with a single peak near the minimum biomass flux.

nbins = 20;
[yUn, xUn] = hist(X2_un(ibm, :), nbins,'linewidth', 2);
[yLims, xLims] = hist(X2_lim(ibm, :), nbins,'linewidth', 2);

f3 = figure;
p1 = plot(xUn, yUn, xLims, yLims,'linewidth', 2);
legend('Aerobic', 'Anaerobic')
xlabel('Biomass Flux (mmol/gDW/h)')
ylabel('# samples')
%% 
% Adding the FVA results to the plot shows that the sampling distributions give 
% more detailed information about the differences between the two models. In particular, 
% we see that the biomass flux minima and maxima are not equally probable. The 
% number of samples from both the aerobic and anaerobic models peaks at the minimum 
% flux of zero, and decreases approximately monotonically towards the maximum 
% growth rate. It decreases more slowly in the unanaerobic model, indicating that 
% higher ATP production is more probable under unlimited oxygen uptake conditions. 
% It is interesting to see that maximum ATP production is highly improbable in 
% both models.

ylims = get(gca, 'ylim');
cUn = get(p1(1), 'color');
cLim = get(p1(2), 'color');

hold on
p2 = plot([minUn(ibm), minUn(ibm)], ylims, 'REPLACE_WITH_DASH_DASH', [maxUn(ibm), maxUn(ibm)], ylims, 'REPLACE_WITH_DASH_DASH','linewidth', 2);
set(p2,'color', cUn)
p3 = plot([minLim(ibm), minLim(ibm)], ylims, 'REPLACE_WITH_DASH_DASH', [maxLim(ibm), maxLim(ibm)], ylims, 'REPLACE_WITH_DASH_DASH','linewidth', 2);
set(p3, 'color', cLim)
legend('Aerobic', 'Anaerobic', 'AerobicFVA','AerobicFVA', 'AnaerobicFVA','AnaerobicFVA')
hold off
%% 
% 
%% Acknowledgements
% Based on a genome-scale sampling tutorial by Hulda S. Haraldsdóttir and German 
% A. Preciat Gonzalez.
%% References
% 1. Orth, J. D., Thiele I., and Palsson, B. Ø. What is flux balance analysis? 
% _Nat. Biotechnol._ 28(3), 245-248 (2010).
% 
% 2. Haraldsdóttir, H. S., Cousins, B., Thiele, I., Fleming, R.M.T., and Vempala, 
% S. CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based 
% metabolic models. _Bioinformatics_. 33(11), 1741-1743 (2016).
% 
% 3. Zhang, Y. and Gao, L. On Numerical Solution of the Maximum Volume Ellipsoid 
% Problem. _SIAM J. Optimiz_. 14(1), 53-76 (2001).
% 
% 4. Berbee, H. C. P., Boender, C. G. E., Rinnooy Ran, A. H. G., Scheffer, C. 
% L., Smith, R. L., Telgen, J. Hit-and-run algorithms for the identification of 
% nonredundant linear inequalities. _Math. Programming_, 37(2), 184-207 (1987).
##### SOURCE END #####
-->
</div></body></html>